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ABSTRACT 

This  paper  is  focused  on  the  dynamic  formulation  of  mechanical  joints  using  different 
approaches  that  lead  to  different  models  with  different  numbers  of  degrees  of  freedom.  Some  of 
these  formulations  allow  for  capturing  the  joint  deformations  using  discrete  elastic  model  while 
the  others  are  continuum-based  and  capture  joint  deformation  modes  that  cannot  be  captured 
using  the  discrete  elastic  joint  models.  Specifically,  three  types  of  joint  formulations  are 
considered  in  this  investigation;  the  ideal ,  compliant  discrete  element,  and  compliant  continuum- 
based  joint  models.  The  ideal  joint  formulation,  which  does  not  allow  for  deformation  degrees  of 
freedom  in  the  case  of  rigid  body  or  small  deformation  analysis,  requires  introducing  a  set  of 
algebraic  constraint  equations  that  can  be  handled  in  computational  multibody  system  (MBS) 
algorithms  using  two  fundamentally  different  approaches:  constrained  dynamics  approach  and 
penalty  method.  When  the  constrained  dynamics  approach  is  used,  the  constraint  equations  must 
be  satisfied  at  the  position,  velocity,  and  acceleration  levels.  The  penalty  method,  on  the  other 
hand,  ensures  that  the  algebraic  equations  are  satisfied  at  the  position  level  only.  In  the  compliant 
discrete  element  joint  formulation,  no  constraint  conditions  are  used;  instead  the  connectivity 
conditions  between  bodies  are  enforced  using  forces  that  can  be  defined  in  their  most  general 
form  in  MBS  algorithms  using  bushing  elements  that  allow  for  the  definition  of  general 
nonlinear  forces  and  moments.  The  new  compliant  continuum-based  joint  formulation,  which  is 
based  on  the  finite  element  (FE)  absolute  nodal  coordinate  formulation  (ANCF),  has  several 
advantages:  (1)  It  captures  modes  of  joint  deformations  that  cannot  be  captured  using  the 
compliant  discrete  joint  models;  (2)  It  leads  to  linear  connectivity  conditions,  thereby  allowing 
for  the  elimination  of  the  dependent  variables  at  a  preprocessing  stage;  (3)  It  leads  to  a  constant 
inertia  matrix  in  the  case  of  chain  like  structure;  and  (4)  It  automatically  captures  the 
deformation  of  the  bodies  using  distributed  inertia  and  elasticity.  The  formulations  of  these  three 
different  joint  models  are  compared  in  order  to  shed  light  on  the  fundamental  differences 
between  them.  Numerical  results  of  a  detailed  tracked  vehicle  model  are  presented  in  order  to 
demonstrate  the  implementation  of  some  of  the  formulations  discussed  in  this  investigation. 

Keywords:  Multibody  system  dynamics;  joint  formulations;  joint  compliance;  penalty  method; 
bushing  elements;  tracked  vehicles;  ANCF  finite  elements. 
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1.  PROBLEM  DEFINITION 

Accurate  formulation  of  mechanical  joints  is  necessary  in  the  computer  simulation  of  multibody 
system  (MBS)  that  represent  many  technological  and  industrial  applications.  An  example  of 
these  MBS  applications,  in  which  accurate  modeling  of  joint  compliance  is  necessary,  is  the 
tracked  vehicle  shown  in  Fig.  1.  The  links  of  the  track  chains  of  this  vehicle  are  connected  by  pin 
joints  that  can  be  subjected  to  significant  stresses  during  the  vehicle  functional  operations. 
Nonetheless,  there  are  different  joint  formulations  that  can  lead  to  different  dynamic  models 
which  have  different  numbers  of  degrees  of  freedom.  This  study  investigates  the  use  of  three 
different  methods  for  formulating  mechanical  joints  in  MBS  applications.  These  three  methods 
are  the  ideal,  the  compliant  discrete  element,  and  the  compliant  continuum-based  joint 
formulations.  These  three  different  methods  are  described  below. 

1 . 1  Ideal  J oint  F ormulation 

The  ideal  joint  formulation  is  based  on  a  set  of  algebraic  equations  that  do  not  account  for  the 
joint  flexibility;  this  is  regardless  of  whether  or  not  the  body  is  flexible.  The  algebraic  joint 
equations  are  expressed  in  terms  of  the  coordinates  of  the  two  bodies  connected  by  the  joint. 
These  algebraic  equations  are  considered  as  constraint  equations  which  can  be  enforced  using 
two  fundamentally  different  methods;  the  constrained  dynamics  approach  and  the  penalty 
method.  In  the  constrained  dynamics  approach,  the  technique  of  Lagrange  multipliers  or  a 
recursive  method  is  used.  In  this  case,  the  joint  constraint  equations  must  be  satisfied  at  the 
position,  velocity,  and  acceleration  levels.  The  number  of  degrees  of  freedom  of  the  model  in  this 
case  is  equal  to  the  number  of  the  system  coordinates  minus  the  number  of  the  algebraic  joint 
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constraint  equations.  In  the  penalty  method,  on  the  other  hand,  the  number  of  degrees  of  freedom 
of  the  model  is  not  affected  by  the  number  of  joint  constraint  equations.  These  joint  algebraic 
equations  are  enforced  using  high  stiffness  penalty  coefficients  that  ensure  that  the  algebraic 
constraint  equations  are  satisfied  at  the  position  level.  The  penalty  method  does  not  ensure  that 
the  constraint  equations  are  satisfied  at  the  velocity  and  acceleration  levels. 

1.2  Compliant  Discrete  Element  Joint  Formulation 

In  this  approach,  no  algebraic  equations  are  used  to  describe  the  joints  between  bodies  in  the 
system.  The  connectivity  between  bodies  is  described  using  force  elements  that  have  forms 
defined  by  the  user  of  the  MBS  code.  MBS  system  codes  have  bushing  elements  that  can  be  used 
to  define  general  linear  or  nonlinear  force  and  moment  expressions.  The  stiffness  and  damping 
coefficients  in  the  force  and  moment  expressions  can  be  selected  by  the  user.  The  busing 
elements  can  be  used  to  model  the  joint  compliance  in  the  case  of  rigid  and  flexible  body 
dynamics.  It  is  important,  however,  to  point  out  that  adding  bushing  elements  has  no  effect  on 
the  number  of  degrees  of  freedom  of  the  model.  Unlike  the  penalty  method,  the  use  of  bushing 
element  does  not  require  the  formulation  of  algebraic  joint  equations.  Bushing  elements  allow  for 
systematically  introducing  three  force  components  and  three  moment  components. 

1.3  Compliant  Continuum-Based  Joint  Formulation 

The  new  finite  element  (FE)  absolute  nodal  coordinate  formulation  (ANCF)  allows  for 
systematically  developing  new  joint  formulations  that  capture  modes  of  deformation  that  cannot 
be  captured  using  the  discrete  joint  models.  It  also  allows  for  modeling  body  flexibility  using 
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new  FE  meshes  that  have  constant  inertia  and  linear  connectivity  conditions.  Specifically,  the 
compliant  ANCF  continuum-based  joint  formulation  has  the  following  advantages: 

1.  ANCF  finite  elements  allow  for  developing  new  joint  formulations  that  capture 
deformation  modes  that  cannot  be  captured  using  compliant  discrete  joint  formulations. 
The  use  of  the  ANCF  gradient  coordinates  allows  for  developing  different  joint  models 
with  different  numbers  of  degrees  of  freedom  that  allow  for  different  strain  modes. 

2.  The  use  of  the  ANCF  gradient  coordinates  allows  for  developing  linear  joint  constraint 
equations.  These  linear  algebraic  equations  can  be  used  to  eliminate  dependent  variables 
at  a  preprocessing  stage,  thereby  significantly  reducing  the  model  dimensionality. 

3.  ANCF  finite  elements  can  also  be  used  to  model  the  body  deformation  in  addition  to  the 
joint  compliance.  Distributed  inertia  and  elasticity  are  used  for  both  body  flexibility  and 
joint  compliance. 

4.  ANCF  finite  elements  lead  to  new  types  of  FE  meshes  that  have  constant  inertia,  a  feature 
that  can  be  exploited  to  develop  a  sparse  matrix  structure  of  the  MBS  dynamic  equations. 

It  is  the  objective  of  this  investigation  to  provide  a  comprehensive  study  of  different  joint 

formulations  and  demonstrate  the  fundamental  differences  between  them  when  applied  to  the 

analysis  of  complex  tracked  vehicle  system  models.  Simulation  models  of  a  typical  tracked 

vehicle  will  be  used  to  compare  the  numerical  results  obtained  using  the  joint  formulations 

discussed  in  this  paper.  These  results  are  obtained  using  the  general  purpose  MBS  computer  code 

SAMS/2000  (Shabana,  2010)  which  allows  for  systematically  modeling  MBS  applications  using 

the  augmented  formulation,  penalty  method,  bushing  elements,  and  ANCF  finite  elements.  In  the 
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augmented  formulation,  the  technique  of  Lagrange  multipliers  is  used.  The  computational 
algorithm  used  in  SAMS/2000  ensures  that  the  algebraic  constraint  equations  are  satisfied  at  the 
position,  velocity,  and  acceleration  levels. 

2.  TRACKED  VEHICLES:  BACKGROUND 

High  mobility  tracked  vehicle  such  as  military  battle  tanks  and  armored  personal  carriers  are 
designed  for  the  mobility  over  rough  and  off-road  terrains.  Investigations  on  the  dynamic 
analysis  of  such  tracked  vehicles  shown  in  Fig.  1  have  been  limited  because  of  the  complexity  of 
the  forces  resulting  from  interaction  between  the  vehicle  components.  These  forces  are  impulsive 
in  nature,  and  their  dynamic  modeling  requires  sophisticated  computational  capabilities.  Several 
two  dimensional  models  for  the  analysis  of  tracked  vehicle  have  been  developed.  Galaitsis 
(1984)  demonstrated  that  the  dynamic  track  tension  and  suspension  loads  in  high  speed  tracked 
vehicles  developed  by  analytical  methods  are  useful  in  evaluating  the  dynamic  characteristics  of 
the  tracked  vehicle  components.  Bando  et  al.  (1991)  outlined  a  procedure  for  the  design  and 
analysis  of  rubber  tracked  small-size  bulldozers,  and  presented  a  computer  simulation  which  was 
used  in  the  evaluation  of  the  vehicle  performance.  Both  steel  and  continuous  rubber  tracks  are 
modeled  by  discretizing  them  into  several  rigid  bodies  connected  by  compliant  elements.  The 
simulation  results  indicate  that  the  vehicle  has  favorable  characteristics,  such  as  less  damage  to 
the  road  surface,  and  reduced  vibration  and  noise.  Murray  and  Canfield  (1992)  used  general 
purpose  multibody  computer  codes  to  model  a  simple  track  link  and  sprocket  system.  The 
behavior  of  the  interaction  between  the  track  link  and  the  sprocket  was  illustrated  graphically  and 
it  was  found  that  the  computer  time  can  be  significantly  reduced  by  using  supercomputers. 
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Nakanishi  et  al.  (1994)  developed  a  two  dimensional  contact  force  model  for  planar  analysis  of 
multibody  tracked  vehicle  systems.  Modal  parameters  (modal  mass,  modal  stiffness,  modal 
damping,  and  mode  shapes),  which  are  determined  experimentally,  are  employed  to  simulate  the 
nonlinear  dynamic  behavior  of  a  multibody  tracked  vehicle  which  consists  of  interconnected 
rigid  and  flexible  components.  The  equations  of  motion  of  the  vehicle  are  formulated  in  terms  of 
a  set  of  modal  and  reference  generalized  coordinates,  and  the  theoretical  basis  for  extracting  the 
component  modal  parameters  of  the  chassis  from  the  modal  parameters  of  the  assembled  vehicle 
is  described. 

A  number  of  approaches  have  been  proposed  in  the  literature  for  developing  three- 
dimensional  MBS  models.  Choi  et  al.  (1998)  developed  the  nonlinear  dynamic  equations  of 
motion  of  the  three-dimensional  multibody  tracked  vehicle  systems,  taking  into  consideration  the 
degrees  of  freedom  of  the  track  chains.  To  avoid  the  solution  of  a  system  of  differential  and 
algebraic  equations,  the  recursive  kinematic  equations  of  the  vehicle  are  expressed  in  terms  of 
the  independent  joint  coordinates.  In  order  to  take  advantage  of  sparse  matrix  algorithms,  the 
independent  differential  equations  of  the  three-dimensional  tracked  vehicles  are  obtained  using 
the  velocity  transformation  method.  Three-dimensional  nonlinear  contact  force  models  that 
describe  the  interaction  between  the  track  links  and  the  vehicle  components  such  as  the  road 
wheels,  sprockets,  and  idlers  as  well  as  the  interaction  between  the  track  links  and  the  ground  are 
developed  and  used  to  define  the  generalized  contact  forces  associated  with  the  vehicle 
generalized  coordinates.  A  computer  simulation  of  a  tracked  vehicle  in  which  the  track  is 
assumed  to  consist  of  track  links  connected  by  a  single  degree  of  freedom  revolute  joint  is 
presented  in  order  to  demonstrate  the  use  of  the  formulations  presented  in  their  study.  Ryu  et  al. 
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(2000)  developed  compliant  track  link  models  and  investigated  the  use  of  these  models  in  the 
dynamic  analysis  of  high-speed,  high-mobility  tracked  vehicles.  The  characteristics  of  the 
compliant  elements  used  in  this  investigation  to  describe  the  track  joints  are  measured 
experimentally.  A  numerical  integration  method  having  a  relatively  large  stability  region  is 
employed  in  order  to  maintain  the  solution  accuracy,  and  a  variable  step  size  integration 
algorithm  is  used  in  order  to  improve  the  efficiency.  The  dimensionality  problem  is  solved  by 
decoupling  the  equations  of  motion  of  the  chassis  and  track  subsystems.  Recursive  methods  are 
used  to  obtain  a  minimum  set  of  equations  for  the  chassis  subsystem.  Several  simulations 
scenarios  including  an  accelerated  motion,  high-speed  motion,  braking,  and  turning  motion  of 
the  high-mobility  vehicle  are  tested  in  order  to  demonstrate  the  effectiveness  and  validity  of  the 
methods  proposed.  Ozaki  and  Shabana  (2003)  evaluated  the  performance  of  different 
formulations  using  a  tracked  vehicle  model  that  is  subjected  to  impulsive  forces.  They  developed 
models  for  joint  constraints  and  the  impulsive  contact  forces  that  result  from  the  interaction 
between  the  track  chains  and  the  vehicle  components  as  well  as  the  interaction  between  these 
chains  and  the  ground.  The  nonlinear  contact  force  models  used  in  their  numerical  study  were 
developed,  and  the  formulations  of  the  generalized  forces  associated  with  the  generalized 
coordinates  used  in  each  of  the  formulations  were  presented.  Ryu  et  al  (2003)  investigated  the 
nonlinear  dynamic  modeling  methods  for  the  virtual  design  of  tracked  vehicles  by  using  MBS 
dynamic  simulation  techniques.  The  results  include  high  oscillatory  signals  resulting  from  the 
impulsive  contact  forces  and  the  use  of  stiff  compliant  elements  to  represent  the  joints  between 
the  track  links.  Each  track  link  is  modeled  as  a  body  which  has  six  degrees  of  freedom,  and  two 
compliant  bushing  elements  are  used  to  connect  track  links.  Efficient  contact  search  and 
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kinematics  algorithms  in  the  context  of  the  compliance  contact  model  are  developed  to  detect  the 
interactions  between  track  links,  road  wheels,  sprockets,  and  ground  for  the  sake  of  speedy  and 
robust  solutions.  Rubinstein  and  Hitron  (2004)  developed  a  three-dimensional  multi-body 
simulation  model  for  simulating  the  dynamic  behavior  of  tracked  off-road  vehicles  using  the 
LMS-DADS  simulation  program.  Each  track  link  is  considered  a  rigid  body  and  is  connected  to 
its  neighboring  track  link  via  a  revolute  joint.  The  road-wheel  track-link  interaction  is  described 
using  three-dimensional  contact  force  elements,  and  the  track-link  terrain  interaction  is  modeled 
using  a  pressure-sinkage  relationship. 

3.  SCOPE  AND  OBJECTIVES  OF  THIS  INVESTIGATION 

While  some  of  the  formulations  used  in  the  analysis  of  tracked  vehicles  require  the  use  of  MBS 
algorithms  that  are  designed  for  solving  systems  of  differential/algebraic  equations  arising  from 
kinematic  joint  constraints,  other  formulations  do  not  require  the  use  of  such  algorithms  but  use 
penalty  forces  or  complaint  elements  to  model  the  chain  dynamics.  The  track  links  are  subjected 
to  high  contact  forces  that  produce  high  stress  levels,  which  often  cause  damage  to  the  track  links 
during  the  functional  operation  of  the  vehicle.  The  objective  of  this  paper  is  to  present  and 
compare  between  different  joint  formulations  that  can  be  used  in  the  dynamic  modeling  of 
complex  tracked  vehicle  systems.  Better  understanding  of  these  formulations  can  lead  to  more 
accurate,  and  possibly  faster,  computer  simulations  that  can  be  the  basis  for  more  reliable 
performance  evaluation  of  the  vehicles.  The  tracked  vehicles  considered  in  this  investigation  are 
assumed  to  consist  of  interconnected  bodies  that  can  have  arbitrary  displacements.  In  the  first 
chain  formulation,  referred  to  in  this  paper  as  the  ideal  joint  chain  model ,  kinematic  joints 
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between  the  track  links  are  described  using  nonlinear  constraint  equations  that  lead  to  significant 
reduction  in  the  number  of  vehicle  degrees  of  freedom.  This  joint  model  does  not  require 
assuming  stiffness  and  damping  for  the  track  link  connectivity,  and  therefore,  does  not  allow  for 
flexibility  between  the  tack  links;  it  requires,  however,  the  solution  of  a  system  of  differential 
and  algebraic  equations  if  redundant  coordinate  formulations  are  used.  Redundant  coordinate 
algorithms  based  on  the  Lagrangian  augmented  form  of  the  equations  of  motion  require  the  use 
of  Newton-Raphson  method  in  order  to  ensure  that  the  constraint  equations  are  satisfied  at  the 
position  level.  Recursive  and  joint  variables  methods  can  also  be  used  instead  of  redundant 
coordinate  formulations  in  order  to  avoid  Newton-Raphson  algorithm.  Another  approach  that  can 
be  used  to  enforce  the  constraint  equations  at  the  position  level  is  the  penalty  method.  This  model 
does  not  lead  to  reduction  in  the  number  of  the  system  degrees  of  freedom. 

Two  other  approaches  that  capture  the  joint  compliance  will  also  be  considered  in  this 
study.  The  first  is  the  compliant  discrete  element  method  that  employs  MBS  bushing  elements  to 
define  the  connectivity  between  the  track  links.  This  approach  as  in  the  case  of  the  penalty 
method  requires  assuming  stiffness  and  damping  coefficients  at  the  connection,  and  therefore,  it 
allows  for  the  flexibility  between  the  track  links.  In  the  second,  the  compliant  continuum-based 
joint  formulation  that  employs  ANCF  finite  elements  is  used.  This  approach,  which  captures  new 
joint  deformation  modes,  leads  to  linear  connectivity  conditions  which  can  be  applied  at  a 
preprocessing  stage  allowing  for  an  efficient  elimination  of  the  dependent  variables,  this  leads  to 
a  constant  inertia  matrix  and  zero  Coriolis  and  centrifugal  forces  (Shabana  et  al.,  2012).  This 
approach  leads  to  new  types  of  FE  chain  meshes  that  have  desirable  characteristics. 
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4.  ALGEBRAIC  CONSTRAINT  EQUATIONS 

In  the  methods  of  constrained  dynamics,  there  are  two  approaches  that  are  often  used  to  model 
ideal  mechanical  joints  that  do  not  account  for  the  effect  of  elasticity  and  damping.  These  two 
methods  are  the  augmented  formulation  that  employs  the  technique  of  Lagrange  multipliers  or 
the  recursive  formulation  which  allows  for  systematic  elimination  of  the  dependent  variables 
using  the  algebraic  equations.  These  two  formulations  are  briefly  discussed  in  this  section. 

4.1  Augmented  Formulation 

In  the  augmented  formulation,  the  constraint  forces  explicitly  appear  in  the  dynamic  equations 
which  are  expressed  in  terms  of  redundant  coordinates.  The  constraint  relationships  are  used  with 
the  differential  equations  of  motion  to  solve  for  the  unknown  accelerations  and  constraint  forces. 
While  this  approach  leads  to  a  sparse  matrix  structure,  it  has  the  drawback  of  increasing  the 
problem  dimensionality  and  it  requires  more  sophisticated  numerical  algorithms  to  solve  the 
resulting  system  of  differential  and  algebraic  equations  (DAE).  Using  the  generalized 
coordinates,  the  equations  of  motion  of  a  body  i  can  be  written  as  (Roberson  and  Schwertassek, 
1988;  Shabana  et  al.,  2008) 

M'q'=Q:+Q:.+Q'  (1) 

. -ij 

where  M'is  the  mass  matrix  of  the  body,  q'  =  R'7  0'7  is  the  vector  of  the  accelerations  of  the 
body  with  R'  defining  the  body  translation  and  0'  defining  the  body  orientation,  Q'  is  the 
vector  of  external  forces,  Q'  is  the  vector  of  the  constraint  forces  which  can  be  written  in  terms 
of  Lagrange  multipliers  X  as  Q'  =  -C7,  X  ,  C  ,  is  the  constraint  Jacobian  matrix  associated  with 
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the  coordinates  of  body  i ,  and  Q';  is  the  vector  of  the  inertia  forces  that  absorb  terms  that  are 
quadratic  in  the  velocities.  The  constraint  equations  at  the  acceleration  level  can  be  written  as 
Cq,q'  =  Q', ,  where  Q';  is  a  vector  that  absorbs  first  derivatives  of  the  coordinates.  Using  Eq.  1 

with  the  constraint  equations  at  the  acceleration  level,  one  obtains 


M  cnrql  to-fQ,/ 
Cq  L 


(2) 


The  matrices  and  vectors  that  appear  in  this  equation  are  the  system  matrices  and  vectors  that  are 
obtained  by  assembling  the  body  matrices  and  vectors.  The  preceding  matrix  equation,  which 
ensures  that  the  constraint  equations  are  satisfied  at  the  acceleration  level,  can  be  solved  for  the 
accelerations  and  Lagrange  multipliers.  In  order  to  ensure  that  the  algebraic  kinematic  constraint 
equations  are  satisfied  at  the  position  and  velocity  levels,  the  independent  accelerations  q,  are 

identified  and  integrated  forward  in  time  in  order  to  determine  the  independent  velocities  q,  and 
independent  coordinates  q  .  Knowing  the  independent  coordinates  from  the  numerical 
integration,  the  dependent  coordinates  qd  can  be  determined  from  the  nonlinear  constraint 
equations  using  an  iterative  Newton-Raphson  algorithm  that  requires  the  solution  of  the  system 
CqAqd  --C,  where  Aq  „  is  the  vector  of  Newton  differences,  and  Cq  is  the  constraint 
Jacobian  matrix  associated  with  the  dependent  coordinates.  Knowing  the  system  coordinates  and 
the  independent  velocities,  the  dependent  velocities  qrf  can  be  determined  by  solving  a  linear 
system  of  algebraic  equations  that  represents  the  constraint  equations  at  the  velocity  level.  This 
linear  system  of  equations  in  the  velocities  can  be  written  as  Cq  qd  =  -Cq  q  -Ct ;  where  Cq  is 
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the  constraint  Jacobian  matrix  associated  with  the  independent  coordinates,  and  Cr  =  dC/dt  is 
the  partial  derivative  of  the  constraint  functions  with  respect  to  time. 

Lagrange  multipliers,  on  the  other  hand,  can  be  used  to  determine  the  constraint  forces. 
For  a  given  joints  ,  the  generalized  constraint  forces  acting  on  body  i ,  connected  by  this  joint, 
can  be  obtained  from  the  equation 

(<>•),  =-(C*)>.=  [>f  tTJ  (3) 

Where  F(  and  T,'  are  the  generalized  joint  forces  associated,  respectively,  with  the  translation 

and  orientation  coordinates  of  body? .  Using  the  results  of  Eq.  3,  the  reaction  forces  at  the  joint 
definition  point  can  be  determined  using  the  concept  of  the  equipollent  systems  of  forces. 

4.2  Recursive  Formulation 

Another  alternate  approach  for  formulating  the  equations  of  motion  of  constrained  mechanical 
systems  is  the  recursive  method,  wherein  the  equations  of  motion  are  formulated  in  terms  of  the 
joint  degrees  of  freedom.  This  formulation  leads  to  a  minimum  set  of  differential  equations  from 
which  the  workless  constraint  forces  are  automatically  eliminated  (Roberson  and  Schwertassek, 
1988;  Shabana,  2010).  The  numerical  procedure  used  in  solving  these  differential  equations  is 
much  simpler  than  the  procedure  used  in  the  solution  of  the  mixed  system  of  differential  and 
algebraic  equations  resulting  from  the  use  of  the  augmented  formulation.  In  the  recursive 
formulation,  the  equations  of  motion  are  formulated  in  terms  of  joint  degrees  of  freedom.  In  this 
formulation,  the  multibody  system  is  assumed  to  consist  of  subsystems,  as  in  the  case  of  the  track 
chains  shown  in  Fig.  2.  The  absolute  coordinates  and  velocities  of  an  arbitrary  body  ?  in  a 
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subsystem  are  expressed  in  terms  of  the  independent  joint  variables  as  well  as  the  absolute 

coordinates  and  velocities  of  body  j  .  If  body  i  is  connected  to  body  j  through  a  revolute  joint, 

which  is  the  case  in  this  subsystem,  the  relative  rotation  is  the  only  degree  of  freedom 

represented  between  the  bodies.  The  connectivity  between  bodies  i  and  body  j  can  then  be 

described  using  the  kinematic  relationships 

R'  +  A  Up  -  R '  -  A  Up  =  0 
to'  =  CD7  +0)''J 

where  R'  is  the  global  position  vector  of  the  origin  of  body  z;  A'  is  the  transformation  matrix 
that  defines  the  body  orientation  and  can  be  expressed  in  terms  of  Euler  parameters;  uj,  and  uj, 
are  the  local  position  vectors  of  point  P  defined  in  the  coordinate  systems  of  body  z  and  j , 
respectively,  to'  and  to7  are,  respectively,  the  absolute  angular  velocity  vectors  of  bodies  z  and 
j  ,  and  o''j  is  the  angular  velocity  vector  of  body  z  with  respect  to  body  j  which  can  be  defined 
as  to1'7  =  '  with  v7  =  A7v7  where  v7  is  a  unit  vector  along  the  axis  of  rotation  defined  in  the 

i 

coordinate  system  of  body  j  ,  and  (j)'  is  the  angle  of  relative  rotation.  By  differentiating  the  first 
equation  in  Eq.  4  twice  and  the  second  once  with  respect  to  time,  one  obtains 

R'  =  -to'  x(to'  xup)-a'  xup  +  R7  +co7  x(oJ  xuj,)  +  a7  xuj, 
a'  =  aJ  +  +  (to7  x  v7 )^' 

In  this  equation,  ak  is  the  absolute  angular  acceleration  vector  of  body  k  .  Using  the  kinematic 
equations  obtained  in  this  section,  one  can  systematically  eliminate  the  dependent  variables  in 
order  to  obtain  a  number  of  differential  equations  of  motion  equal  to  the  number  of  the  system 
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degrees  of  freedom.  Using  this  approach,  one  obtains  a  dense  inertia  matrix  in  a  system  of 
dynamic  equations  that  does  not  have  constraint  forces.  A  second,  alternative  approach  is  to  use 
the  kinematic  equations  developed  in  this  section  to  determine  all  the  system  absolute 
coordinates  and  velocities.  One  can  then  construct  Eq.  2,  which  can  be  solved  for  the 
accelerations  and  Lagrange  multipliers.  Using  the  absolute  acceleration  relationships  of  Eq.  5, 
one  can  determine  the  relative  joint  accelerations.  The  joint  accelerations  can  be  integrated 
forward  in  time  in  order  to  determine  the  joint  coordinates  and  velocities. 


5.  GENERALIZED  FORCES 

In  defining  the  joint  forces  between  the  track  links,  it  is  important  to  understand  the  relationship 
and  differences  between  the  generalized  and  the  Cartesian  moments  (Roberson  and 
Schwertassek,  1988;  Shabana,  2010).  This  is  important  in  interpreting  the  reaction  forces  of  the 
ideal  joints  and  also  important  in  the  implementation  of  the  penalty  method  and  bushing 
elements.  Let  F'  be  a  force  vector  that  acts  at  a  point  P'  on  a  rigid  body  i.  If  this  force  vector  is 
assumed  to  be  defined  in  the  global  coordinate  system,  then  the  virtual  work  of  this  force  vector 

■T 

can  be  written  as  SWle  =  F'  Sr'p ,  where  S r'p  can  be  found  using  the  virtual  change  in  the  position 
vector  of  an  arbitrary  point  on  rigid  body  i  as 


Sr ■;  = 


50 

[i  -A  uj,G  ] 

1 

!  ^ 
_ i 

(6) 


In  this  equation.  A' is  the  transformation  matrix  that  defines  the  body  orientation,  is  the  skew 
symmetric  matrix  associated  with  the  vector  uj,  that  defines  the  local  coordinates  of  the  point 
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P' ,  and  G'  is  the  matrix  that  relates  the  angular  velocity  vector  defined  in  the  body 
coordinate  system  to  the  time  derivatives  of  the  orientation  coordinates,  that  is  =  G  O  .  Note 
that  since  up  =  A'uPA'  ,  Eq.  6  can  be  written  as  Jrp  =  JR'  — u'pG' J0'  Using  this  equation  in 

■T  ■  -T 

the  virtual  work  expression,  one  obtains  SW'  =F'  JR'  -F'  u'pG'J0' ,  which  can  be  written  as 

SWl  =F;  JR'+F''J0'  (7) 

■T  . T 

where  Fp  =  F' ,  and  F'0  =  -G'  upF' .  These  equations  imply  that  a  force  that  acts  at  an  arbitrary 
point  on  the  rigid  body  i  is  equipollent  to  another  system  defined  at  the  reference  point  that 

■T  T 

consists  of  the  same  force  and  a  set  of  generalized  forces,  defined  by  F(  =  -G'  upF'  associated 
with  the  orientation  coordinates  of  the  body  (Roberson  and  Schwertassek,  1988;  Shabana,  2010). 

•  -T 

Since  up  is  a  skew-symmetric  matrix,  it  follows  that  u'p  =  -ulp  .  Using  this  identity,  one  can 

■T  ... 

show  that  the  generalized  moment  can  be  written  as  Fp  =G'  (u),  x  F' )  ,  where  M'(  =upxF'  is 

the  Cartesian  moment  resulting  from  the  application  of  the  force  F',  and  G'  is  the  matrix  that 
relates  the  angular  velocity  vector  to'  defined  in  the  global  coordinate  system  to  the  time 
derivatives  of  the  orientation  coordinates,  that  is  to'  =  G'0' .  It  follows  that  the  relationship 

■T 

between  the  generalized  and  Cartesian  moment  is  Fp  =G'  M'.  If  the  components  of  the 
moment  vector  are  defined  in  the  body  coordinate  system,  one  has  F'd=G'  M'( ,  where 
M';  =  A'  M'(.  The  relationships  developed  in  this  section  will  be  used  in  the  formulation  of 
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the  joint  forces  in  the  case  of  the  penalty  method.  These  relationships  will  also  be  used  in  the 
computer  implementation  of  the  bushing  element  in  general  MBS  algorithms. 

6.  JOINT  FORMULATIONS 

In  this  paper,  four  different  joint  models  are  considered;  two  models  are  based  on  algebraic 
equations  that  require  the  use  of  the  methods  of  constrained  dynamics  or  the  penalty  method.  In 
the  case  of  constrained  dynamics,  an  alternative  to  the  use  of  Lagrange  multipliers  is  the  use  of 
the  recursive  methods,  as  previously  discussed.  When  Lagrange  multiplier  technique  or  the 
recursive  methods  are  used,  the  constraint  equations  must  be  satisfied  at  the  position,  velocity, 
and  acceleration  levels.  The  penalty  method,  on  the  other  hand,  satisfies  the  algebraic  constraint 
equations  at  the  position  level  only. 

6.1  Constraint  Equations 

The  revolute  (pin)  joint  is  used  in  this  section  as  an  example  to  demonstrate  the  formulation  of 
the  algebraic  constraint  equations.  This  joint  has  been  used  in  the  literature  in  the  modeling  of 
the  track  chains.  As  shown  in  Fig.  2,  the  track  chain  can  be  assumed  to  consist  of  links  connected 
to  each  other  by  a  pin  joint  that  allows  for  the  relative  rotation  between  them.  In  general  MBS 
algorithms,  the  nonlinear  algebraic  equations  that  define  the  pin  joint  are  expressed  in  terms  of 
the  absolute  coordinates  of  the  two  bodies  i  and  j  connected  by  the  joint.  The  five  algebraic 
constraint  equations  that  eliminate  five  degrees  of  freedom  can  be  written  in  terms  of  the 
absolute  Cartesian  coordinates  of  the  two  bodies  as  (Shabana,  2010) 
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(8) 


where  v'  and  v7  are  two  vectors  defined  along  the  joint  axis  on  bodies  i  and  j ,  respectively; 


v' ,  v',Vj  form  an  orthogonal  triad  defined  on  body  i ;  k0  is  a  constant  (Shabana,  2010);  and 


4  =r‘p-  rj,  =  R  +  A'iip  -Ry  -  AJuJP  (9) 

In  this  equation,  A'  and  Aj  are  the  transformation  matrices  that  define  the  orientation  of  bodies 
i  and  j ,  respectively,  and  uj,  and  uj,  are  the  local  position  vectors  of  points  P'  and  P1  with 
respect  to  bodies  i  and  j  ,  respectively.  Points  P‘  and  Pj  are  defined  on  the  axis  of  the  pin  joint 
on  bodies  i  and  j ,  respectively.  One  can  show  that  the  Jacobian  matrix  of  the  pin  joint 
constraints  is  defined  as 


v;Th; 

v'fH7 

vtH2 

v'l'H' 

h;  +  v'7h'p 

-v'fH ' 

h;wx 

-v'TH  Jp 

2rfH'p 

-2rfH  Jp 

(10) 


where  C  ,  and  C  ,  are  the  constraint  Jacobian  matrices  associated  with  the  coordinates  of  bodies 
q  q 

i  and  j ,  respectively;  and  other  vectors  and  matrices  that  appear  in  the  preceding  equation  are 


»  1  A'b'/G'],HJ=[i  A^G'],H;=|f  =  A(A'v;).H'2=^  =  ^-(A'r2), 


dw\ 


Hy,=  — =— (AV). 

dqJ  8qJ  V  ; 
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Another  alternate  approach  for  formulating  the  revolute  joint  constraints  is  to  consider  it  as  a 
special  case  of  the  spherical  joint  in  which  the  relative  rotation  between  the  two  bodies  is 
allowed  only  along  the  joint  axis.  If  point  P  is  the  joint  definition  point,  and  v'  and  \J  are  two 
vectors  defined  along  the  joint  axis  on  bodies  i  and  j  ,  respectively,  the  constraint  equations  of 

the  revolute  joint  can  be  written  as  C(q',q;  )  =  r ,JP  v'fv;  v^v7]  =  0 .  The  last  two  equations 

in  this  equation  guarantee  that  the  two  vectors  v'  and  \J  remain  parallel,  thereby  eliminating  the 
freedom  of  the  relative  rotation  between  the  two  bodies  in  two  perpendicular  directions. 

6.2  Penalty  Method 

The  constraint  equations  that  describe  the  connectivity  between  the  track  links  can  be  enforced 
using  the  penalty  approach.  In  this  case,  these  algebraic  equations  are  satisfied  only  at  the 
position  level.  The  penalty  method  does  not  lead  to  elimination  of  degrees  of  freedom,  and 
therefore,  it  is  conceptually  different  from  the  case  of  Lagrange  multiplier  technique  or  the 
recursive  approach.  In  order  to  demonstrate  the  penalty  approach,  the  violations  in  the  constraint 
equations  of  a  revolute  joint  k  can  be  written  as 

dr=[rp  (11) 

Using  this  violation  dk ,  a  restoring  force  vector  can  then  be  defined  as  fk  =  kk  dk  +ckdk  ,  where 
kk ,  and  ck  are  assumed  penalty  stiffness  and  damping  coefficients,  respectively,  and  dt  is  the 
time  derivative  of  the  violation  vector  d,  .  The  virtual  work  of  this  restoring  force  f;.  can  then  be 
written  as  5Wk  =  -fkSdk ,  which  can  be  written  as 
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SW?  =-^-^(y^J)-V(v^) 


(12) 


where  Srl,p  =  SR'  - u^G  ^0  -  SRj  +  u' G  W ,  FB  =  kkriJp  +  ck if ,  S(\[T\j )  =  \jT S\[  +  v[T S\j , 


S(y'2t\j)  =  \jT S\‘2  +  y'2tSyj  ,  Fl  =kk  (v'fv')  +  c; 


cl  (v'JV) 


dt 


,  and  F2  =  kk(\iT2\j)  +  cl 


d(  yiT2yj) 


dt 


Equation  12  can  be  used  to  define  a  set  of  generalized  forces  acting  on  bodies  i  and  j  that 
maintain  the  connectivity  between  the  two  bodies  as 


SW’’  =  Q‘RrSR'  +  Q'lSQ'  +  Q''SR'+QJ7Sd! 


(13) 


Which  can  be  rewritten  in  a  compact  form  as  5W  =  Q'/b'q'  +  Q'/  SqJ ,  where  q'  and  q'  are  the 
generalized  coordinates  of  bodies  i  and  j  ,  and 
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(14) 


where  M'  =  ,  M '  =  Fjvj  v; ,  M',  =  -Fjv'2  v' ,  and  M(  =  FjV(  \j .  Equation  14  defines  the 

generalized  forces  associated  with  the  absolute  Cartesian  coordinates  due  to  the  revolute  joint 
connection  between  bodies  i  and  j .  With  a  proper  selection  of  the  penalty  coefficients,  these 
forces  ensure  that  the  constraint  equations  are  satisfied  at  the  position  level. 


7.  COMPLIANT  DISCRETE  ELEMENT  JOINT  FORMULATION 
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The  compliant  discrete  element  joint  formulation  allows  for  introducing  joint  deformations.  In 
this  approach,  no  algebraic  equations  are  enforced.  In  most  MBS  computer  codes,  the  compliant 
discrete  element  joint  method  can  be  applied  using  the  standard  MBS  bushing  element  that 
allows  for  introducing  three  force  and  three  moment  components  that  can  be  linear  or  nonlinear 
functions  of  the  body  coordinates.  As  shown  in  Fig.  3,  the  position  vectors  np  and  u Jp  of  two 

points,  FJ7  and  Pj  on  body  j,  can  be  used  to  define  one  axis  of  the  coordinate  system  of  the 
bushing  element  as  n7  =  (il^  —  u  ^  j  j  |u^  —  | ,  where  n7  is  one  of  the  bushing  axes  defined  in 

the  body  j  coordinate  system.  This  axis  can  then  be  defined  in  the  global  coordinate  system  as 

n7  =  A7n7  ,  where  A7  is  the  transformation  matrix  that  defines  the  orientation  of  the  coordinate 
system  of  body  j  in  the  global  system.  Using  this  axis,  one  can  define  the  directional  properties 
of  the  bushing  element  with  the  other  two  axes  of  the  bushing  coordinate  system  being  defined 
using  the  transformation  matrix  A1’1  =  Tt/  t/  n7J  where  t(  and  V2  are  the  two  unit  vectors 

that  complete  the  three  orthogonal  axes  of  the  bushing  element  coordinate  system.  Assuming  that 
body  j  is  a  rigid  body,  the  bushing  coordinate  system  can  be  defined  with  respect  to  the  global 

coordinate  system  as  Ahj  =  AjAbj .  Choosing  points  P'  and  P/  to  initially  coincide;  one  can 
define  the  bushing  deformation  and  rate  of  deformation  vectors  in  the  bushing  coordinate  system 
as  5hij  =  Ahj  r'J ,  and  6hij  =  Ahj  r'7,  respectively,  where  r'7  =  r',  -  r7  is  the  position  vector  of 
point  Pxj  with  respect  to  point  P‘ . 

The  rotational  deformation  of  the  bushing  element  can  be  obtained  using  the  transformation 
matrix  that  defines  the  orientation  of  the  bushing  coordinate  system  on  body  i  with  respect  to  the 
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bushing  coordinate  system  on  body  j  .  This  matrix  is  defined  as  AbiJ  =  Ahj  Abl ,  where  A hj  is  the 
orientation  matrix  of  the  bushing  coordinate  system  on  body  j ,  while  Ahl  is  the  orientation 
matrix  of  the  bushing  coordinate  system  with  respect  to  the  coordinate  system  of  body  i  that  is 
defined  as  Abl  =A'Abl.  Assuming  that  the  relative  rotations  between  bodies  /  and  j  arc  small, 
the  relative  rotation  matrix,  AbiJ  can  be  used  to  extract  three  relative  rotations  defined  in  the 
bushing  coordinate  system,  Qbij  —  \j)  bij  0  blJ  ().  hT^  .  The  relative  angular  velocity  between  the 

two  bodies  defined  in  the  bushing  coordinate  system  can  also  be  written  as  (Obij  =  Ah]  (co'  —  to7) , 
where  to'  and  to7  are  the  absolute  angular  velocity  vectors  of  bodies  i  and  j ,  respectively, 
defined  in  the  global  coordinate  system.  The  bushing  stiffness  and  damping  coefficients  are  often 
determined  using  experimental  testing,  and  these  coefficients  are  defined  generally  in  the 
bushing  coordinate  system.  Let  K;.  and  C,  be  the  translational  stiffness  and  damping  matrices, 
respectively,  defined  with  respect  to  the  bushing  coordinate  system;  and  assume  that  the 
rotational  stiffness  and  damping  matrices  are  K  0  and  C.g,  respectively.  In  terms  of  translational 

and  rotational  stiffness  and  damping  matrices,  the  force  vector  defined  in  the  bushing  coordinate 
system  can  be  written  as 


Fb]_\Kr  0  Jdbu~\  \Cr  0Tg^ 

M*J“Lo  kJHIo  cJ[©*b 


(15) 


This  force  vector  can  then  be  defined  in  the  global  coordinate  system,  and  the  results  can  be  used 
to  define  the  generalized  bushing  forces  and  moments  acting  on  the  two  bodies  as  previously 
described  in  this  paper. 


21 


UNCLASSIFIED 


8.  COMPLIANT  CONTINUUM-BASED  JOINT  FORMULATION 

The  compliant  continuum-based  joint  formulation  allows  capturing  joint  strain  modes  that  cannot 
be  captured  using  the  compliant  discrete  element  joint  method.  When  ANCF  finite  elements  are 
used,  one  can  develop  new  FE  meshes  that  have  linear  connectivity  and  constant  inertia 
(Shabana  et  al.,  2012).  This  allows  for  systematically  eliminating  dependent  variables  at  a 
preprocessing  stage,  and  as  a  result,  there  is  no  need  for  the  use  of  joint  formulations  in  the  main 
processor.  This  approach  can  be  used  to  develop  new  spatial  chain  models  where  the  modes  of 
deformations  at  the  definition  points  of  the  joints  that  allow  for  rigid  body  rotations  between 
ANCF  finite  elements  can  be  captured.  The  displacement  field  of  an  ANCF  finite  element,  as  the 
one  shown  in  Fig.  4,  can  be  written  as  r(x,  y,z,t)  =  s(x, y,z)e(t)  where  x,y,  and  z  are  the 
element  spatial  coordinates;  t  is  time;  S  is  the  element  shape  function  matrix,  and  e  is  the 
vector  of  element  nodal  coordinates.  Using  this  displacement  field,  the  equations  of  a  pin  joint 
between  elements  i  and  j  can  be  written  using  the  six  scalar  equations  r'  =  rj ,  r'a  =  r j  ,  where 
a  is  the  coordinate  line  that  defines  the  joint  axis;  a  can  be  x,  y ,  or  z  or  any  other  coordinate 
line.  The  six  scalar  equations  eliminate  six  degrees  of  freedom;  three  translations,  two  rotations, 
and  one  deformation  mode.  Therefore,  this  joint  has  five  modes  of  deformation  that  include 
stretch  and  shear  modes.  This  ANCF  revolute  joint  model  ensures  C1  continuity  with  respect  to 
the  coordinate  line  a  and  C°  continuity  with  respect  to  the  other  two  parameters.  It  follows  that 
the  Lagrangian  strain  component  saa  =  (Viy  -l)/2  is  continuous  at  the  joint  definition  point, 

while  the  other  five  strain  components  can  be  discontinuous.  The  resulting  joint  constraint 

equations  are  linear,  and  therefore,  can  be  applied  at  a  preprocessing  stage  to  systematically 

eliminate  the  dependent  variables.  Using  these  equations,  one  can  develop  a  new  kinematically 
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linear  FE  mesh  for  flexible-link  chains  in  which  the  links  can  have  arbitrarily  large  relative 
rotations. 

In  the  numerical  investigation  presented  in  this  paper,  a  three-dimensional  cable  element 
is  used  to  model  the  flexibility  of  the  chain  links.  For  the  three-dimensional  cable  element  used 
in  the  compliant  continuum-based  model,  the  nodal  coordinates  can  be  written  as 


e7  = 


r7  (x(  =0)  r7  (x(  =  0)  r7  (x(  =  /)  r7  (x(  =  /)]  ,  with  the  element  shape  function 


defined  as  S7  =  [.s’, I  .s2I  .s,I  .s4l]  where  I  is  the  2x2  identity  matrix  (Gerstmayr  and 
Shabana,  2006).  The  shape  functions  5,  ,  for  i  =  1,  2,  3,  4,  are  defined  as 


s^l-3^  +  2^,  s2=l(4- 2f~+Z3) 


s3= 


(16) 


where  E,  -  x(  / /  and  /  is  the  length  of  the  element.  The  constraint  forces,  QJC  ,  can  then  be  found 
using  the  equation  of  motion  for  a  single  body  (element)  j  as  M7e7  =(Q7  +Q;(,„ +Q;!)  +  Qc> 
where  Q7 ,  Q70„ ,  and  Q7  are  the  elastic  forces,  forces  of  contacts  between  the  chain  links  and 
other  bodies  in  the  system,  and  gravity  forces,  respectively,  associated  with  the  element  nodal 
coordinates  e7  and  M7  is  the  symmetric  mass  matrix  of  the  finite  element  j  defined  as 

M  =  |  pjSjTSj  dVj ,  where  p '  is  the  mass  density  of  the  material  points  of  the  element  and  V  J 

vJ 

is  the  element  volume.  Since  the  shape  function  for  a  three-dimensional  cable  element  only 
depends  on  the  xl  spatial  coordinate,  the  mass  matrix  can  then  be  defined  as 
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d  h 
~2  2  l 


MJ  -  |  1 1 Po&j  S Jdxldx2dx3  -d'h'  p^S'  S' dxl 


-d  -h  0 


(17) 


For  the  cable  element,  this  mass  matrix  can  be  written  explicitly  as 


M  .//!/• 
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where  di  ,hj ,  and  are  the  width,  height,  and  initial  density,  respectively,  of  element  j.  It  is 

important  to  note  that  the  mass  matrix  is  constant  in  both  two-dimensional  and  three-dimensional 
cases  and  leads  to  zero  centrifugal  and  Coriolis  forces  when  the  body  experiences  an  arbitrary 
large  deformation  and  finite  rotation.  The  virtual  work  of  the  elastic  forces  can  then  be  defined  as 


/  / 

SWs  =  J  EAs]  ,&!  dx\  +|  EIkSkcIx^ 

0  0 


(19) 


where  E  is  the  modulus  of  elasticity,  A  is  the  element  cross  section,  I  is  the  second  moment  of 
area  and  k  is  the  curvature.  The  elastic  forces  of  the  cable  element  can  also  be  evaluated  using 
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i  i 

the  following  expression  for  the  strain  energy  U  =  J  EA(sn)2dx1  +J  EI(ic)2dxx  as  Q  s 


dU_ 

V  3e  j 


Similarly,  the  gravity  forces  can  be  found  using  the  virtual  work  of  the  gravity  forces, 
We  =  J[0  0  -  pQg^StdV .  Substituting  in  the  shape  functions  and  integrating  over  the 


volume  leads  to  SW  =  -mg 


00-00  —  00-00  — 

2  12  2  12 


St  ,  which  leads  to 


Q,  =  -mg 


00-00  —  00-00  — 

2  12  2  12 


(20) 


Using  the  principle  of  virtual  work  with  the  mass  matrix  and  nodal  accelerations  as  previously 
defined,  the  joint  forces  for  a  specific  element  can  then  be  easily  calculated. 

9.  SIMULATION  RESULTS 


In  this  section,  numerical  results,  obtained  using  the  tracked  vehicle  model  shown  in  Fig.  1,  are 
used  to  compare  the  different  joint  formulations  presented  in  this  investigation.  The  tracked 
vehicle  modeled  is  an  armored  personnel  carrier  that  consists  of  a  chassis,  idler,  sprocket,  5  road- 
wheels,  and  64  track  links  on  each  track  side  (right  and  left).  Figure  5  shows  the  engagement  of 
the  track  links  with  some  of  the  vehicle  components,  while  more  details  about  the  road  wheel 
arrangement  and  the  configuration  of  the  suspension  system  of  the  tracked  vehicle  model  are 
shown  in  Fig.  6.  Table  1  shows  the  inertia  properties  for  all  the  tracked  vehicle  model 
components  used  in  this  simulation.  More  specifications  of  this  vehicle  can  be  obtained  from 
open  sources  (Ml  13,  2003).  The  vehicle  has  a  suspension  system  that  consists  of  road  arms 
placed  between  the  road  wheels  and  chassis  as  well  as  shock  absorbers  connected  to  each  road 
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arm.  Table  2  shows  the  stiffness  and  the  damping  coefficients  of  the  contact  models  used.  In  this 
study,  two  different  simulation  scenarios,  one  with  the  suspension  system  and  the  other  without 
the  suspension  system,  will  be  considered  to  study  the  effect  of  using  the  suspension  system  on 
the  results.  The  road  arms  and  the  sprockets  are  connected  to  the  chassis  by  revolute  joints,  and 
the  road  arms  are  connected  to  the  road-wheels  by  revolute  joints.  The  track  links  are  connected 
to  each  other  using  revolute  joints,  which  can  be  modeled  using  the  constraint  equations,  penalty 
method,  bushing  element,  or  ANCF  finite  elements  as  previously  mentioned.  Tensioners  are 
added  to  the  system  by  connecting  each  idler  to  a  tensioner  with  a  revolute  joint  and  connecting 
the  tensioner  to  the  chassis  with  a  prismatic  joint  to  ensure  only  translation  between  them.  The 
angular  velocities  of  the  sprockets  of  the  vehicle  model  considered  in  this  numerical 
investigation  are  assumed  to  increase  after  1  sec  until  they  reach  a  constant  value  of  25  rad/sec 
after  8  seconds  as  shown  in  Fig.  7.  Figure  8  shows  the  chassis  vertical  displacement  using  the 
joint  model  with  and  without  the  suspension  system.  The  results  presented  in  this  figure  show 
that  the  model  with  a  suspension  system  allows  for  more  vibration  (a  7  cm  initial  drop  compared 
to  0.64  cm),  which  helps  make  the  model  more  realistic  as  compared  to  the  previously 
unsuspended  model.  The  ANCF  finite  element  for  the  track  links  used  in  this  investigation  is  a 
three-dimensional  steel  cable  element  with  a  modulus  of  rigidity  of  76.9  GPa,  a  Young’s 
modulus  of  200  GPa,  and  a  mass  density  of  7.80x10  kg/m  .  Figures  9  and  10  show, 
respectively,  the  chassis  forward  position  and  velocity  results  obtained  using  different  joint 
models.  While  the  results  presented  in  these  figures  show  good  agreement,  the  computational 
time  varies  when  these  different  models  are  used.  The  constrained  joint  model  takes  less 
computational  time  compared  to  the  penalty,  bushing  element,  and  ANCF  joint  models;  the 
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penalty  model  CPU  time  is  six  times  that  of  the  constraint  model  CPU  time,  while  the  bushing 
model  CPU  time  is  three  times,  and  the  ANCF  joint  model  CPU  time  is  5  times  that  of  the 
constraint  model  CPU  time.  This  increase  in  the  penalty  and  bushing  element  CPU  time  is 
attributed  to  the  high  stiffness  coefficients  used  in  both  of  these  models.  Efforts  are  currently 
being  made  to  improve  the  efficiency  of  the  ANCF  model.  Figure  11  shows  the  motion  trajectory 
of  a  track  link  in  the  chassis  coordinate  system  using  the  constraint  model,  the  penalty  model,  the 
bushing  element  model,  and  the  ANCF  joint  model,  respectively.  The  results  in  this  figure  show 
good,  although  not  exact,  agreement  between  the  different  models.  Figures  12-15  show  the  joint 
forces  in  the  longitudinal  and  the  vertical  directions,  respectively,  obtained  using  the  constrained, 
penalty,  and  ANCF  joint  models.  Very  high  frequencies  were  filtered  out  in  all  models  using  a 
low  pass  FFT  with  a  cut-off  frequency  of  30  Hz  in  order  to  show  clearly  the  nominal  values  of 
the  joint  forces  presented  in  Figs.  12b  and  13b.  The  results  presented  in  these  figures  are 
obtained  using  a  stiffness  coefficient  of  lxlO9  N/m  and  damping  coefficient  of  lxlO5  N.s/m  for 
the  penalty  method  model.  The  results  show  good  agreement  between  the  constrained  and 
penalty  joint  models  with  the  ANCF  joint  model,  as  expected,  showing  lower  force  magnitude 
due  to  the  flexibility  of  the  elements.  Figures  14  and  15  show  the  same  results  in  the  case  of  a 
stiffness  coefficient  of  lxlO7  N/m  and  damping  coefficient  of  lxlO5  N.s/m  for  the  penalty 
method  model.  The  results  of  Figs.  14  and  15,  which  show  significant  differences  between  the 
two  ideal  joint  models,  demonstrate  the  drawback  of  the  penalty  method  when  the  penalty 
stiffness  coefficient  is  reduced.  Similar  results  can  be  expected  in  the  case  of  the  bushing  element 
models  where  the  forces  obtained  also  depend  on  the  stiffness  and  damping  coefficient  of  the 
joint.  Figure  16  shows  the  joint  deformation  predicted  using  the  penalty  model  using  different 
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stiffness  coefficients.  The  penalty  model  with  a  stiffness  coefficient  of  lxlO9  N/m  shows  much 
less  deformation,  less  than  0.075  mm,  while  the  penalty  model  with  a  stiffness  of  lxlO7  N/m  has 
much  more  deformation,  over  1.5  mm,  between  the  track  links. 

10.  SUMMARY  AND  CONCULSIONS 

In  this  investigation,  different  MBS  joint  formulations  are  presented  and  compared  using  detailed 
tracked  vehicle  models.  Three  main  joint  formulations  are  discussed;  they  are  the  ideal  joint 
formulation,  the  compliant  discrete  element  joint  formulation,  and  the  compliant  continuum- 
based  joint  formulation.  The  ideal  joint  formulation  is  developed  to  eliminate  the  relative 
displacement  between  the  two  bodies  connected  by  the  joint.  This  can  be  achieved  by  enforcing  a 
set  of  joint  algebraic  equations  using  a  constrained  dynamics  approach  or  by  using  the  penalty 
method.  The  constrained  dynamics  approach  eliminates  degrees  of  freedom  and  ensures  that  the 
constraint  equations  are  satisfied  at  the  position,  velocity,  and  acceleration  levels.  The  penalty 
method,  on  the  other  hand,  does  not  reduce  the  number  of  degrees  of  freedom  and  ensures  that 
the  constraint  equations  are  satisfied  at  the  position  level  only  provided  that  a  high  stiffness 
coefficient  is  used.  The  compliant  discrete  element  formulation,  which  allows  for  joint 
deformations,  can  be  systematically  applied  using  a  standard  MBS  bushing  element  that  allows 
for  six  degrees  of  freedom  of  relative  motion.  The  compliant  continuum-based  approach  can  be 
used  to  develop  new  joints  that  capture  deformation  modes  that  are  not  captured  by  the 
compliant  discrete  element  joint  formulation.  ANCF  finite  elements  can  be  used  to 
systematically  develop  new  joints  with  distributed  elasticity  and  linear  connectivity  conditions. 


28 


UNCLASSIFIED 


As  discussed  in  this  paper,  it  is  important  to  choose  the  proper  stiffness  and  damping 
coefficients  when  the  penalty  method  and  the  bushing  elements  are  used.  Numerical  results 
were  presented  in  order  to  compare  between  different  methods.  The  ideal  joint  formulation 
produces  the  desired  joint  kinematics  and  accurate  joint  forces.  The  same  is  true  with  the 
penalty  force  based  joint  when  large  penalty  stiffness  coefficients  are  used.  Penalty  force 
based  joint  construction  has  been  shown  to  be  sensitive  to  the  selection  of  penalty  stiffness  with 
the  higher  stiffness  coefficients  leading  to  better  overall  results.  However,  higher  penalty 
stiffness  increases  CPU  time  significantly  due  to  higher  frequencies.  The  penalty  method  and 
bushing  element  models  each  have  much  larger  CPU  times  than  the  ideal  constrained  model  due 
to  these  high  stiffness  coefficients.  The  results  presented  in  this  investigation  show  that  the 
ANCF  joint  model  leads  to  lower  force  predictions  which  can  be  attributed  to  the  track  link 
flexibility. 
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Table  1:  Mass  and  inertia  values  for  tracked  vehicle  parts 


Part 

Mass  (kg) 

Ixx  ( kg.m2) 

Iyy  ( kg.m2) 

Izz  ( kg.m2) 

Ixy,  Iyz,  Ixz  (kg.m”) 

Chassis 

5489.2 

1786.9 

10450 

10721 

0 

Sprocket 

436.67 

13.868 

12.216 

12.216 

0 

Idler 

429.57 

14.698 

12.545 

12.545 

0 

Road  Wheel 

561.07 

26.063 

19.819 

19.819 

0 

Road  Arm 

75.264 

0.77085 

0.37632 

0.77085 

0 

Track  Link 

18.024 

0.04113 

0.22463 

0.25256 

0 

Table  2:  Contact  Parameters 


Parameters 

Sprocket-Track 

Contact 

Road  Wheel-Track 

Contact 

Ground-Track 

Contact 

k 

2.00xl0b  N/m 

2.00xl0&  N/m 

2.00xl0&  N/m 

c 

5.00xl0J  N-s/m 

5-OOxlO3  N-s/m 

5.00xl03  N-s/m 

F 

0.150 

0.100 

0.300 
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Chassis 


Road- wheel 


Figure  1:  The  SAMS/2000  tracked  vehicle  model 


Figure  2:  Revolute  joint 
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Figure  3:  Bushing  element 


Figure  4:  ANCF  beam  element  coordinates 
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Figure  5:  Tracked  vehicle  contact  models 


Track  tensioner 


Shock  absorber  Road  ann 


Figure  6:  Suspension  system  layout  of  the  tracked  vehicle 
(www.combatreform.org/SoucyBandTrackInstallationA-654RB.pdt~) 
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Figure  7 :  Sprocket  angular  velocity 


Figure  8:  Chassis  vertical  displacement 
(  A  with  suspension,  — ▼ —  without  suspension) 
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Figure  9:  Chassis  forward  position 

(  A  Constrained  joint  model,  — ▼ —  Penalty  method  model,  ■  Bushing  element  model, 

•  ANCF  model) 


Figure  10:  Chassis  forward  velocity 

(  A  Constrained  joint  model,  — ▼ —  Penalty  method  model,  ■  Bushing  element  model, 

•  ANCF  model) 
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Figure  11:  Trajectory  motion  of  a  track  link  in  the  chassis  coordinate  system 
(  A  Constrained  joint  model,  — ▼ —  Penalty  method  model,  ■  Bushing  element  model, 

— ANCF  model) 


8 


40000 
30000  - 
20000- 
10000  - 

-10000- 
-20000  - 
-30000  - 
-40000 


-10000- 
-20000  - 
-30000  - 
-40000 


Time  (s) 


Figure  12a:  Joint  longitudinal  forces 

(— A-  Constrained  joint  model,  — Penalty  method  model  k  =  1  x  109  N/m,  — ANCF  model) 
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Figure  12b:  Filtered  joint  longitudinal  forces  using  FFT 
(— A-  Constrained  joint  model,  — Penalty  method  model  k  =  1  x  109  N/m,  — ANCF  model) 


Figure  13a:  Joint  vertical  forces 

(— *—  Constrained  joint  model,  — Penalty  method  model  k  =  lxlO9  N/m.  — ANCF  model) 
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Figure  13b:  Filtered  joint  vertical  forces  using  FFT 
(— *—  Constrained  joint  model,  — Penalty  method  model  k  =  lxlO9  N/m.  — ANCF  model) 


Figure  14:  Joint  Longitudinal  forces 

Constrained  joint  model,  — Penalty  method  model  k  =  lxlO7  N/m,  —■—ANCF  model) 
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Figure  15:  Joint  vertical  forces 

(— A-  Constrained  joint  model,  — Penalty  method  model  k  =  lxlO7  N/m,  — ANCF  model) 


Figure  16:  Joint  deformation  using  penalty  model 
(-A-  k  =  1  x  109  N/m,  k  =  lx  107  N/m) 
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